home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / hamradio / sgp4_pl2.zip / SGP4SDP4.PAS < prev    next >
Pascal/Delphi Source File  |  1992-09-03  |  41KB  |  1,036 lines
  1. Unit SGP4SDP4;
  2. {           Author:  Dr TS Kelso }
  3. { Original Version:  1991 Oct 30}
  4. { Current Revision:  1992 Sep 03}
  5. {          Version:  1.50 }
  6. {        Copyright:  1991-1992, All Rights Reserved }
  7. {$N+}
  8.  
  9. INTERFACE
  10.   Uses Support,
  11.        SGP_Init,
  12.        SGP_Math,SGP_Time;
  13.  
  14. Procedure SGP(time : double;
  15.        var pos,vel : vector);
  16. Procedure SGP4(tsince : double;
  17.             var iflag : integer;
  18.           var pos,vel : vector);
  19. Procedure SDP4(tsince : double;
  20.             var iflag : integer;
  21.           var pos,vel : vector);
  22.  
  23. IMPLEMENTATION
  24.   Uses SGP_Intf;
  25.  
  26. var
  27. {dpinit}
  28.   eqsq,siniq,cosiq,rteqsq,ao,cosq2,sinomo,cosomo,
  29.   bsq,xlldot,omgdt,xnodot,xnodp                   : double;
  30. {dpsec/dpper}
  31.   xll,omgasm,xnodes,_em,xinc,xn,t                 : double;
  32.   qoms2t                                          : double;
  33.  
  34. Procedure Define_Derived_Constants;
  35.   begin
  36.   xke    := Sqrt(3600*ge/Cube(xkmper));  {Sqrt(ge) ER^3/min^2}
  37.   qoms2t := Sqr(Sqr(qo-s));              {(qo-s)^4 ER^4}
  38.   end; {Define_Derived_Constants}
  39.  
  40. Procedure SGP4(tsince : double;
  41.             var iflag : integer;
  42.           var pos,vel : vector);
  43.   label
  44.     9,10,90,100,110,130,140;
  45.   const
  46.     a1     : double = 0;  a3ovk2 : double = 0;  ao     : double = 0;
  47.     aodp   : double = 0;  aycof  : double = 0;  betao  : double = 0;
  48.     betao2 : double = 0;  c1     : double = 0;  c1sq   : double = 0;
  49.     c2     : double = 0;  c3     : double = 0;  c4     : double = 0;
  50.     c5     : double = 0;  coef   : double = 0;  coef1  : double = 0;
  51.     cosio  : double = 0;  d2     : double = 0;  d3     : double = 0;
  52.     d4     : double = 0;  del1   : double = 0;  delmo  : double = 0;
  53.     delo   : double = 0;  eeta   : double = 0;  eosq   : double = 0;
  54.     eta    : double = 0;  etasq  : double = 0;  isimp  : integer = 0;
  55.     omgcof : double = 0;  omgdot : double = 0;  perige : double = 0;
  56.     pinvsq : double = 0;  psisq  : double = 0;  qoms24 : double = 0;
  57.     s4     : double = 0;  sinio  : double = 0;  sinmo  : double = 0;
  58.     t2cof  : double = 0;  t3cof  : double = 0;  t4cof  : double = 0;
  59.     t5cof  : double = 0;  temp   : double = 0;  temp1  : double = 0;
  60.     temp2  : double = 0;  temp3  : double = 0;  theta2 : double = 0;
  61.     theta4 : double = 0;  tsi    : double = 0;  x1m5th : double = 0;
  62.     x1mth2 : double = 0;  x3thm1 : double = 0;  x7thm1 : double = 0;
  63.     xhdot1 : double = 0;  xlcof  : double = 0;  xmcof  : double = 0;
  64.     xmdot  : double = 0;  xnodcf : double = 0;  xnodot : double = 0;
  65.     xnodp  : double = 0;
  66.   var
  67.     i                                                  : integer;
  68.     cosuk,sinuk,rfdotk,vx,vy,vz,ux,uy,uz,xmy,xmx,
  69.     cosnok,sinnok,cosik,sinik,rdotk,xinck,xnodek,uk,rk,
  70.     cos2u,sin2u,u,sinu,cosu,betal,rfdot,rdot,r,pl,elsq,
  71.     esine,ecose,epw,temp6,temp5,temp4,cosepw,sinepw,
  72.     capu,ayn,xlt,aynl,xll,axn,xn,beta,xl,e,a,tfour,
  73.     tcube,delm,delomg,templ,tempe,tempa,xnode,tsq,xmp,
  74.     omega,xnoddf,omgadf,xmdf,x,y,z,xdot,ydot,zdot      : double;
  75.   begin
  76. { Recover original mean motion (xnodp) and semimajor axis (aodp) }
  77. { from input elements. }
  78.   if (iflag = 0) then
  79.     goto 100;
  80.   a1 := Power(xke/xno,tothrd);
  81.   cosio := Cos(xincl);
  82.   theta2 := cosio*cosio;
  83.   x3thm1 := 3*theta2 - 1;
  84.   eosq := eo*eo;
  85.   betao2 := 1 - eosq;
  86.   betao := sqrt(betao2);
  87.   del1 := 1.5*ck2*x3thm1/(a1*a1*betao*betao2);
  88.   ao := a1*(1 - del1*(0.5*tothrd + del1*(1 + 134/81*del1)));
  89.   delo := 1.5*ck2*x3thm1/(ao*ao*betao*betao2);
  90.   xnodp := xno/(1 + delo);
  91.   aodp := ao/(1 - delo);
  92. { Initialization }
  93. { For perigee less than 220 kilometers, the isimp flag is set and
  94.   the equations are truncated to linear variation in sqrt a and
  95.   quadratic variation in mean anomaly.  Also, the c3 term, the
  96.   delta omega term, and the delta m term are dropped. }
  97.   isimp := 0;
  98.   if (aodp*(1 - eo)/ae) < (220/xkmper + ae) then
  99.     isimp := 1;
  100. { For perigee below 156 km, the values of s and qoms2t are altered. }
  101.   s4 := s;
  102.   qoms24 := qoms2t;
  103.   perige := (aodp*(1 - eo) - ae)*xkmper;
  104.   if (perige >= 156) then
  105.     goto 10;
  106.   s4 := perige - 78;
  107.   if (perige > 98) then
  108.     goto 9;
  109.   s4 := 20;
  110.   9:
  111.   qoms24 := Power((120 - s4)*ae/xkmper,4);
  112.   s4 := s4/xkmper + ae;
  113.  10:
  114.   pinvsq := 1/(aodp*aodp*betao2*betao2);
  115.   tsi := 1/(aodp - s4);
  116.   eta := aodp*eo*tsi;
  117.   etasq := eta*eta;
  118.   eeta := eo*eta;
  119.   psisq := abs(1 - etasq);
  120.   coef := qoms24*Power(tsi,4);
  121.   coef1 := coef/Power(psisq,3.5);
  122.   c2 := coef1*xnodp*(aodp*(1 + 1.5*etasq + eeta*(4 + etasq))
  123.       + 0.75*ck2*tsi/psisq*x3thm1*(8 + 3*etasq*(8 + etasq)));
  124.   c1 := bstar*c2;
  125.   sinio := Sin(xincl);
  126.   a3ovk2 := -xj3/ck2*Power(ae,3);
  127.   c3 := coef*tsi*a3ovk2*xnodp*ae*sinio/eo;
  128.   x1mth2 := 1 - theta2;
  129.   c4 := 2*xnodp*coef1*aodp*betao2*(eta*(2 + 0.5*etasq)
  130.       + eo*(0.5 + 2*etasq) - 2*ck2*tsi/(aodp*psisq)
  131.       *(-3*x3thm1*(1 - 2*eeta + etasq*(1.5 - 0.5*eeta))
  132.       + 0.75*x1mth2*(2*etasq - eeta*(1 + etasq))*Cos(2*omegao)));
  133.   c5 := 2*coef1*aodp*betao2*(1 + 2.75*(etasq + eeta) + eeta*etasq);
  134.   theta4 := theta2*theta2;
  135.   temp1 := 3*ck2*pinvsq*xnodp;
  136.   temp2 := temp1*ck2*pinvsq;
  137.   temp3 := 1.25*ck4*pinvsq*pinvsq*xnodp;
  138.   xmdot := xnodp + 0.5*temp1*betao*x3thm1
  139.          + 0.0625*temp2*betao*(13 - 78*theta2 + 137*theta4);
  140.   x1m5th := 1 - 5*theta2;
  141.   omgdot := -0.5*temp1*x1m5th + 0.0625*temp2*(7 - 114*theta2 +395*theta4)
  142.           + temp3*(3 - 36*theta2 + 49*theta4);
  143.   xhdot1 := -temp1*cosio;
  144.   xnodot := xhdot1 + (0.5*temp2*(4 - 19*theta2)
  145.           + 2*temp3*(3 - 7*theta2))*cosio;
  146.   omgcof := bstar*c3*Cos(omegao);
  147.   xmcof := -tothrd*coef*bstar*ae/eeta;
  148.   xnodcf := 3.5*betao2*xhdot1*c1;
  149.   t2cof := 1.5*c1;
  150.   xlcof := 0.125*a3ovk2*sinio*(3 + 5*cosio)/(1 + cosio);
  151.   aycof := 0.25*a3ovk2*sinio;
  152.   delmo := Power(1 + eta*Cos(xmo),3);
  153.   sinmo := Sin(xmo);
  154.   x7thm1 := 7*theta2 - 1;
  155.   if (isimp = 1) then
  156.     goto 90;
  157.   c1sq := c1*c1;
  158.   d2 := 4*aodp*tsi*c1sq;
  159.   temp := d2*tsi*c1/3;
  160.   d3 := (17*aodp + s4)*temp;
  161.   d4 := 0.5*temp*aodp*tsi*(221*aodp + 31*s4)*c1;
  162.   t3cof := d2 + 2*c1sq;
  163.   t4cof := 0.25*(3*d3 + c1*(12*d2 + 10*c1sq));
  164.   t5cof := 0.2*(3*d4 + 12*c1*d3 + 6*d2*d2 + 15*c1sq*(2*d2 + c1sq));
  165.  90:
  166.   iflag := 0;
  167. { Update for secular gravity and atmospheric drag. }
  168. 100:
  169.   xmdf := xmo + xmdot*tsince;
  170.   omgadf := omegao + omgdot*tsince;
  171.   xnoddf := xnodeo + xnodot*tsince;
  172.   omega := omgadf;
  173.   xmp := xmdf;
  174.   tsq := tsince*tsince;
  175.   xnode := xnoddf + xnodcf*tsq;
  176.   tempa := 1 - c1*tsince;
  177.   tempe := bstar*c4*tsince;
  178.   templ := t2cof*tsq;
  179.   if (isimp = 1) then
  180.     goto 110;
  181.   delomg := omgcof*tsince;
  182.   delm := xmcof*(Power(1 + eta*Cos(xmdf),3) - delmo);
  183.   temp := delomg + delm;
  184.   xmp := xmdf + temp;
  185.   omega := omgadf - temp;
  186.   tcube := tsq*tsince;
  187.   tfour := tsince*tcube;
  188.   tempa := tempa - d2*tsq - d3*tcube - d4*tfour;
  189.   tempe := tempe + bstar*c5*(Sin(xmp) - sinmo);
  190.   templ := templ + t3cof*tcube + tfour*(t4cof + tsince*t5cof);
  191. 110:
  192.   a := aodp*Sqr(tempa);
  193.   e := eo - tempe;
  194.   xl := xmp + omega + xnode + xnodp*templ;
  195.   beta := sqrt(1 - e*e);
  196.   xn := xke/Power(a,1.5);
  197. { Long period periodics }
  198.   axn := e*Cos(omega);
  199.   temp := 1/(a*beta*beta);
  200.   xll := temp*xlcof*axn;
  201.   aynl := temp*aycof;
  202.   xlt := xl + xll;
  203.   ayn := e*Sin(omega) + aynl;
  204. { Solve Kepler's Equation }
  205.   capu := fmod2p(xlt - xnode);
  206.   temp2 := capu;
  207.   for i := 1 to 10 do
  208.     begin
  209.     sinepw := Sin(temp2);
  210.     cosepw := Cos(temp2);
  211.     temp3 := axn*sinepw;
  212.     temp4 := ayn*cosepw;
  213.     temp5 := axn*cosepw;
  214.     temp6 := ayn*sinepw;
  215.     epw := (capu - temp4 + temp3 - temp2)/(1 - temp5 - temp6) + temp2;
  216.     if (abs(epw - temp2) <= e6a) then
  217.       goto 140;
  218. 130:
  219.     temp2 := epw;
  220.     end; {for i}
  221. { Short period preliminary quantities }
  222. 140:
  223.   ecose := temp5 + temp6;
  224.   esine := temp3 - temp4;
  225.   elsq := axn*axn + ayn*ayn;
  226.   temp := 1 - elsq;
  227.   pl := a*temp;
  228.   r := a*(1 - ecose);
  229.   temp1 := 1/r;
  230.   rdot := xke*sqrt(a)*esine*temp1;
  231.   rfdot := xke*sqrt(pl)*temp1;
  232.   temp2 := a*temp1;
  233.   betal := sqrt(temp);
  234.   temp3 := 1/(1 + betal);
  235.   cosu := temp2*(cosepw - axn + ayn*esine*temp3);
  236.   sinu := temp2*(sinepw - ayn - axn*esine*temp3);
  237.   u := actan(sinu,cosu);
  238.   sin2u := 2*sinu*cosu;
  239.   cos2u := 2*cosu*cosu - 1;
  240.   temp := 1/pl;
  241.   temp1 := ck2*temp;
  242.   temp2 := temp1*temp;
  243. { Update for short periodics }
  244.   rk := r*(1 - 1.5*temp2*betal*x3thm1) + 0.5*temp1*x1mth2*cos2u;
  245.   uk := u - 0.25*temp2*x7thm1*sin2u;
  246.   xnodek := xnode + 1.5*temp2*cosio*sin2u;
  247.   xinck := xincl + 1.5*temp2*cosio*sinio*cos2u;
  248.   rdotk := rdot - xn*temp1*x1mth2*sin2u;
  249.   rfdotk := rfdot + xn*temp1*(x1mth2*cos2u + 1.5*x3thm1);
  250. { Orientation vectors }
  251.   sinuk := Sin(uk);
  252.   cosuk := Cos(uk);
  253.   sinik := Sin(xinck);
  254.   cosik := Cos(xinck);
  255.   sinnok := Sin(xnodek);
  256.   cosnok := Cos(xnodek);
  257.   xmx := -sinnok*cosik;
  258.   xmy := cosnok*cosik;
  259.   ux := xmx*sinuk + cosnok*cosuk;
  260.   uy := xmy*sinuk + sinnok*cosuk;
  261.   uz := sinik*sinuk;
  262.   vx := xmx*cosuk - cosnok*sinuk;
  263.   vy := xmy*cosuk - sinnok*sinuk;
  264.   vz := sinik*cosuk;
  265. { Position and velocity }
  266.   x := rk*ux;  pos[1] := x;
  267.   y := rk*uy;  pos[2] := y;
  268.   z := rk*uz;  pos[3] := z;
  269.   xdot := rdotk*ux + rfdotk*vx;  vel[1] := xdot;
  270.   ydot := rdotk*uy + rfdotk*vy;  vel[2] := ydot;
  271.   zdot := rdotk*uz + rfdotk*vz;  vel[3] := zdot;
  272.   end; {Procedure SGP4}
  273.  
  274. Procedure Deep(ideep : integer);
  275.   const
  276.     zns    =  1.19459E-5;     c1ss   =  2.9864797E-6;   zes    =  0.01675;
  277.     znl    =  1.5835218E-4;   c1l    =  4.7968065E-7;   zel    =  0.05490;
  278.     zcosis =  0.91744867;     zsinis =  0.39785416;     zsings = -0.98088458;
  279.     zcosgs =  0.1945905;      zcoshs =  1;              zsinhs =  0;
  280.     q22    =  1.7891679E-6;   q31    =  2.1460748E-6;   q33    =  2.2123015E-7;
  281.     g22    =  5.7686396;      g32    =  0.95240898;     g44    =  1.8014998;
  282.     g52    =  1.0508330;      g54    =  4.4108898;      root22 =  1.7891679E-6;
  283.     root32 =  3.7393792E-7;   root44 =  7.3636953E-9;   root52 =  1.1428639E-7;
  284.     root54 =  2.1765803E-9;   thdt   =  4.3752691E-3;
  285.   label
  286.     {dpinit}
  287.     5,10,20,30,40,45,50,55,60,65,70,80,
  288.     {dpsec}
  289.     90,100,105,110,120,125,130,140,150,152,154,160,165,170,175,180,
  290.     {dpper}
  291.     205,210,218,220,230;
  292. const { Typed constants to retain values between ENTRY calls }
  293.     iresfl : integer = 0;  isynfl : integer = 0;
  294.     iret   : integer = 0;  iretn  : integer = 0;
  295.     ls     : integer = 0;
  296.     a1     : double = 0;    a2     : double = 0;    a3     : double = 0;
  297.     a4     : double = 0;    a5     : double = 0;    a6     : double = 0;
  298.     a7     : double = 0;    a8     : double = 0;    a9     : double = 0;
  299.     a10    : double = 0;    ainv2  : double = 0;    alfdp  : double = 0;
  300.     aqnv   : double = 0;    atime  : double = 0;    betdp  : double = 0;
  301.     bfact  : double = 0;    c      : double = 0;    cc     : double = 0;
  302.     cosis  : double = 0;    cosok  : double = 0;    cosq   : double = 0;
  303.     ctem   : double = 0;    d2201  : double = 0;    d2211  : double = 0;
  304.     d3210  : double = 0;    d3222  : double = 0;    d4410  : double = 0;
  305.     d4422  : double = 0;    d5220  : double = 0;    d5232  : double = 0;
  306.     d5421  : double = 0;    d5433  : double = 0;    dalf   : double = 0;
  307.     day    : double = 0;    dbet   : double = 0;    del1   : double = 0;
  308.     del2   : double = 0;    del3   : double = 0;    delt   : double = 0;
  309.     dls    : double = 0;    e3     : double = 0;    ee2    : double = 0;
  310.     eoc    : double = 0;    eq     : double = 0;    f2     : double = 0;
  311.     f220   : double = 0;    f221   : double = 0;    f3     : double = 0;
  312.     f311   : double = 0;    f321   : double = 0;    f322   : double = 0;
  313.     f330   : double = 0;    f441   : double = 0;    f442   : double = 0;
  314.     f522   : double = 0;    f523   : double = 0;    f542   : double = 0;
  315.     f543   : double = 0;    fasx2  : double = 0;    fasx4  : double = 0;
  316.     fasx6  : double = 0;    ft     : double = 0;    g200   : double = 0;
  317.     g201   : double = 0;    g211   : double = 0;    g300   : double = 0;
  318.     g310   : double = 0;    g322   : double = 0;    g410   : double = 0;
  319.     g422   : double = 0;    g520   : double = 0;    g521   : double = 0;
  320.     g532   : double = 0;    g533   : double = 0;    gam    : double = 0;
  321.     omegaq : double = 0;    pe     : double = 0;    pgh    : double = 0;
  322.     ph     : double = 0;    pinc   : double = 0;    pl     : double = 0;
  323.     preep  : double = 0;    s1     : double = 0;    s2     : double = 0;
  324.     s3     : double = 0;    s4     : double = 0;    s5     : double = 0;
  325.     s6     : double = 0;    s7     : double = 0;    savtsn : double = 0;
  326.     se     : double = 0;    se2    : double = 0;    se3    : double = 0;
  327.     sel    : double = 0;    ses    : double = 0;    sgh    : double = 0;
  328.     sgh2   : double = 0;    sgh3   : double = 0;    sgh4   : double = 0;
  329.     sghl   : double = 0;    sghs   : double = 0;    sh     : double = 0;
  330.     sh2    : double = 0;    sh3    : double = 0;    sh1    : double = 0;
  331.     shs    : double = 0;    si     : double = 0;    si2    : double = 0;
  332.     si3    : double = 0;    sil    : double = 0;    sini2  : double = 0;
  333.     sinis  : double = 0;    sinok  : double = 0;    sinq   : double = 0;
  334.     sinzf  : double = 0;    sis    : double = 0;    sl     : double = 0;
  335.     sl2    : double = 0;    sl3    : double = 0;    sl4    : double = 0;
  336.     sll    : double = 0;    sls    : double = 0;    sse    : double = 0;
  337.     ssg    : double = 0;    ssh    : double = 0;    ssi    : double = 0;
  338.     ssl    : double = 0;    stem   : double = 0;    step2  : double = 0;
  339.     stepn  : double = 0;    stepp  : double = 0;    temp   : double = 0;
  340.     temp1  : double = 0;    thgr   : double = 0;    x1     : double = 0;
  341.     x2     : double = 0;    x2li   : double = 0;    x2omi  : double = 0;
  342.     x3     : double = 0;    x4     : double = 0;    x5     : double = 0;
  343.     x6     : double = 0;    x7     : double = 0;    x8     : double = 0;
  344.     xfact  : double = 0;    xgh2   : double = 0;    xgh3   : double = 0;
  345.     xgh4   : double = 0;    xh2    : double = 0;    xh3    : double = 0;
  346.     xi2    : double = 0;    xi3    : double = 0;    xl     : double = 0;
  347.     xl2    : double = 0;    xl3    : double = 0;    xl4    : double = 0;
  348.     xlamo  : double = 0;    xldot  : double = 0;    xli    : double = 0;
  349.     xls    : double = 0;    xmao   : double = 0;    xnddt  : double = 0;
  350.     xndot  : double = 0;    xni    : double = 0;    xno2   : double = 0;
  351.     xnodce : double = 0;    xnoi   : double = 0;    xnq    : double = 0;
  352.     xomi   : double = 0;    xpidot : double = 0;    xqncl  : double = 0;
  353.     z1     : double = 0;    z11    : double = 0;    z12    : double = 0;
  354.     z13    : double = 0;    z2     : double = 0;    z21    : double = 0;
  355.     z22    : double = 0;    z23    : double = 0;    z3     : double = 0;
  356.     z31    : double = 0;    z32    : double = 0;    z33    : double = 0;
  357.     zcosg  : double = 0;    zcosgl : double = 0;    zcosh  : double = 0;
  358.     zcoshl : double = 0;    zcosi  : double = 0;    zcosil : double = 0;
  359.     ze     : double = 0;    zf     : double = 0;    zm     : double = 0;
  360.     zmo    : double = 0;    zmol   : double = 0;    zmos   : double = 0;
  361.     zn     : double = 0;    zsing  : double = 0;    zsingl : double = 0;
  362.     zsinh  : double = 0;    zsinhl : double = 0;    zsini  : double = 0;
  363.     zsinil : double = 0;    zx     : double = 0;    zy     : double = 0;
  364.   begin
  365.   case ideep of
  366.     dpinit : begin { Entrance for deep space initialization }
  367.              thgr := Thetag(epoch);
  368.              eq := eo;
  369.              xnq := xnodp;
  370.              aqnv := 1/ao;
  371.              xqncl := xincl;
  372.              xmao := xmo;
  373.              xpidot := omgdt + xnodot;
  374.              sinq := Sin(xnodeo);
  375.              cosq := Cos(xnodeo);
  376.              omegaq := omegao;
  377.     { Initialize lunar solar terms }
  378.           5: day := ds50 + 18261.5;  {Days since 1900 Jan 0.5}
  379.              if (day = preep) then
  380.                Goto 10;
  381.              preep := day;
  382.              xnodce := 4.5236020 - 9.2422029E-4*day;
  383.              stem := Sin(xnodce);
  384.              ctem := Cos(xnodce);
  385.              zcosil := 0.91375164 - 0.03568096*ctem;
  386.              zsinil := Sqrt(1 - zcosil*zcosil);
  387.              zsinhl := 0.089683511*stem/zsinil;
  388.              zcoshl := Sqrt(1 - zsinhl*zsinhl);
  389.              c := 4.7199672 + 0.22997150*day;
  390.              gam := 5.8351514 + 0.0019443680*day;
  391.              zmol := fmod2p(c - gam);
  392.              zx := 0.39785416*stem/zsinil;
  393.              zy := zcoshl*ctem + 0.91744867*zsinhl*stem;
  394.              zx := Actan(zx,zy);
  395.              zx := gam + zx - xnodce;
  396.              zcosgl := Cos(zx);
  397.              zsingl := Sin(zx);
  398.              zmos := 6.2565837 + 0.017201977*day;
  399.              zmos := fmod2p(zmos);
  400.     { Do solar terms }
  401.          10: savtsn := 1E20;
  402.              zcosg := zcosgs;
  403.              zsing := zsings;
  404.              zcosi := zcosis;
  405.              zsini := zsinis;
  406.              zcosh := cosq;
  407.              zsinh := sinq;
  408.              cc := c1ss;
  409.              zn := zns;
  410.              ze := zes;
  411.              zmo := zmos;
  412.              xnoi := 1/xnq;
  413.              ls := 30; {assign 30 to ls}
  414.          20: a1 := zcosg*zcosh + zsing*zcosi*zsinh;
  415.              a3 := -zsing*zcosh + zcosg*zcosi*zsinh;
  416.              a7 := -zcosg*zsinh + zsing*zcosi*zcosh;
  417.              a8 := zsing*zsini;
  418.              a9 := zsing*zsinh + zcosg*zcosi*zcosh;
  419.              a10 := zcosg*zsini;
  420.              a2 := cosiq*a7 +  siniq*a8;
  421.              a4 := cosiq*a9 +  siniq*a10;
  422.              a5 := -siniq*a7 +  cosiq*a8;
  423.              a6 := -siniq*a9 +  cosiq*a10;
  424.              x1 := a1*cosomo + a2*sinomo;
  425.              x2 := a3*cosomo + a4*sinomo;
  426.              x3 := -a1*sinomo + a2*cosomo;
  427.              x4 := -a3*sinomo + a4*cosomo;
  428.              x5 := a5*sinomo;
  429.              x6 := a6*sinomo;
  430.              x7 := a5*cosomo;
  431.              x8 := a6*cosomo;
  432.              z31 := 12*x1*x1 - 3*x3*x3;
  433.              z32 := 24*x1*x2 - 6*x3*x4;
  434.              z33 := 12*x2*x2 - 3*x4*x4;
  435.              z1 := 3*(a1*a1 + a2*a2) + z31*eqsq;
  436.              z2 := 6*(a1*a3 + a2*a4) + z32*eqsq;
  437.              z3 := 3*(a3*a3 + a4*a4) + z33*eqsq;
  438.              z11 := -6*a1*a5 + eqsq*(-24*x1*x7 - 6*x3*x5);
  439.              z12 := -6*(a1*a6 + a3*a5)
  440.                   + eqsq*(-24*(x2*x7 + x1*x8) - 6*(x3*x6 + x4*x5));
  441.              z13 := -6*a3*a6 + eqsq*(-24*x2*x8 - 6*x4*x6);
  442.              z21 := 6*a2*a5 + eqsq*(24*x1*x5 - 6*x3*x7);
  443.              z22 := 6*(a4*a5 + a2*a6)
  444.                   + eqsq*(24*(x2*x5 + x1*x6) - 6*(x4*x7 + x3*x8));
  445.              z23 := 6*a4*a6 + eqsq*(24*x2*x6 - 6*x4*x8);
  446.              z1 := z1 + z1 + bsq*z31;
  447.              z2 := z2 + z2 + bsq*z32;
  448.              z3 := z3 + z3 + bsq*z33;
  449.              s3 := cc*xnoi;
  450.              s2 := -0.5*s3/rteqsq;
  451.              s4 := s3*rteqsq;
  452.              s1 := -15*eq*s4;
  453.              s5 := x1*x3 + x2*x4;
  454.              s6 := x2*x3 + x1*x4;
  455.              s7 := x2*x4 - x1*x3;
  456.              se := s1*zn*s5;
  457.              si := s2*zn*(z11 + z13);
  458.              sl := -zn*s3*(z1 + z3 - 14 - 6*eqsq);
  459.              sgh := s4*zn*(z31 + z33 - 6);
  460.              sh := -zn*s2*(z21 + z23);
  461.              if (xqncl < 5.2359877E-2) then
  462.                sh := 0;
  463.              ee2 := 2*s1*s6;
  464.              e3 := 2*s1*s7;
  465.              xi2 := 2*s2*z12;
  466.              xi3 := 2*s2*(z13 - z11);
  467.              xl2 := -2*s3*z2;
  468.              xl3 := -2*s3*(z3 - z1);
  469.              xl4 := -2*s3*(-21 - 9*eqsq)*ze;
  470.              xgh2 := 2*s4*z32;
  471.              xgh3 := 2*s4*(z33 - z31);
  472.              xgh4 := -18*s4*ze;
  473.              xh2 := -2*s2*z22;
  474.              xh3 := -2*s2*(z23 - z21);
  475.              case ls of
  476.                30 : Goto 30;
  477.                40 : Goto 40;
  478.              else
  479.                     Halt;
  480.              end; {case}
  481.     { Do lunar terms }
  482.          30: sse := se;
  483.              ssi := si;
  484.              ssl := sl;
  485.              ssh := sh/siniq;
  486.              ssg := sgh - cosiq*ssh;
  487.              se2 := ee2;
  488.              si2 := xi2;
  489.              sl2 := xl2;
  490.              sgh2 := xgh2;
  491.              sh2 := xh2;
  492.              se3 := e3;
  493.              si3 := xi3;
  494.              sl3 := xl3;
  495.              sgh3 := xgh3;
  496.              sh3 := xh3;
  497.              sl4 := xl4;
  498.              sgh4 := xgh4;
  499.              zcosg := zcosgl;
  500.              zsing := zsingl;
  501.              zcosi := zcosil;
  502.              zsini := zsinil;
  503.              zcosh := zcoshl*cosq + zsinhl*sinq;
  504.              zsinh := sinq*zcoshl - cosq*zsinhl;
  505.              zn := znl;
  506.              cc := c1l;
  507.              ze := zel;
  508.              zmo := zmol;
  509.              ls := 40; {assign 40 to ls}
  510.              Goto 20;
  511.          40: sse := sse + se;
  512.              ssi := ssi + si;
  513.              ssl := ssl + sl;
  514.              ssg := ssg + sgh - cosiq/siniq*sh;
  515.              ssh := ssh + sh/siniq;
  516.     { Geopotential resonance initialization for 12 hour orbits }
  517.              iresfl := 0;
  518.              isynfl := 0;
  519.              if (xnq < 0.0052359877) and (xnq > 0.0034906585) then
  520.                Goto 70;
  521.              if (xnq < 8.26E-3) or (xnq > 9.24E-3) then
  522.                exit;
  523.              if (eq < 0.5) then
  524.                exit;
  525.              iresfl := 1;
  526.              eoc := eq*eqsq;
  527.              g201 := -0.306 - (eq - 0.64)*0.440;
  528.              if (eq > 0.65) then
  529.                Goto 45;
  530.              g211 := 3.616 - 13.247*eq + 16.290*eqsq;
  531.              g310 := -19.302 + 117.390*eq - 228.419*eqsq + 156.591*eoc;
  532.              g322 := -18.9068 + 109.7927*eq - 214.6334*eqsq + 146.5816*eoc;
  533.              g410 := -41.122 + 242.694*eq - 471.094*eqsq + 313.953*eoc;
  534.              g422 := -146.407 + 841.880*eq - 1629.014*eqsq + 1083.435*eoc;
  535.              g520 := -532.114 + 3017.977*eq - 5740*eqsq + 3708.276*eoc;
  536.              Goto 55;
  537.          45: g211 := -72.099 + 331.819*eq - 508.738*eqsq + 266.724*eoc;
  538.              g310 := -346.844 + 1582.851*eq - 2415.925*eqsq + 1246.113*eoc;
  539.              g322 := -342.585 + 1554.908*eq - 2366.899*eqsq + 1215.972*eoc;
  540.              g410 := -1052.797 + 4758.686*eq - 7193.992*eqsq + 3651.957*eoc;
  541.              g422 := -3581.69 + 16178.11*eq - 24462.77*eqsq + 12422.52*eoc;
  542.              if (eq > 0.715) then
  543.                Goto 50;
  544.              g520 := 1464.74 - 4664.75*eq + 3763.64*eqsq;
  545.              Goto 55;
  546.          50: g520 := -5149.66 + 29936.92*eq - 54087.36*eqsq + 31324.56*eoc;
  547.          55: if (eq >= (0.7)) then
  548.                Goto 60;
  549.              g533 := -919.2277 + 4988.61*eq - 9064.77*eqsq + 5542.21*eoc;
  550.              g521 := -822.71072 + 4568.6173*eq - 8491.4146*eqsq + 5337.524*eoc;
  551.              g532 := -853.666 + 4690.25*eq - 8624.77*eqsq + 5341.4*eoc;
  552.              Goto 65;
  553.          60: g533 := -37995.78 + 161616.52*eq - 229838.2*eqsq + 109377.94*eoc;
  554.              g521 := -51752.104 + 218913.95*eq - 309468.16*eqsq + 146349.42*eoc;
  555.              g532 := -40023.88 + 170470.89*eq - 242699.48*eqsq + 115605.82*eoc;
  556.          65: sini2 := siniq*siniq;
  557.              f220 := 0.75*(1 + 2*cosiq + cosq2);
  558.              f221 := 1.5*sini2;
  559.              f321 := 1.875*siniq*(1 - 2*cosiq - 3*cosq2);
  560.              f322 := -1.875*siniq*(1 + 2*cosiq - 3*cosq2);
  561.              f441 := 35*sini2*f220;
  562.              f442 := 39.3750*sini2*sini2;
  563.              f522 := 9.84375*siniq*(sini2*(1 - 2*cosiq - 5*cosq2)
  564.                    + 0.33333333*(-2 + 4*cosiq + 6*cosq2));
  565.              f523 := siniq*(4.92187512*sini2*(-2 - 4*cosiq + 10*cosq2)
  566.                    + 6.56250012*(1 + 2*cosiq - 3*cosq2));
  567.              f542 := 29.53125*siniq*(2 - 8*cosiq + cosq2*(-12 + 8*cosiq + 10*cosq2));
  568.              f543 := 29.53125*siniq*(-2 - 8*cosiq + cosq2*(12 + 8*cosiq - 10*cosq2));
  569.              xno2 := xnq*xnq;
  570.              ainv2 := aqnv*aqnv;
  571.              temp1 := 3*xno2*ainv2;
  572.              temp := temp1*root22;
  573.              d2201 := temp*f220*g201;
  574.              d2211 := temp*f221*g211;
  575.              temp1 := temp1*aqnv;
  576.              temp := temp1*root32;
  577.              d3210 := temp*f321*g310;
  578.              d3222 := temp*f322*g322;
  579.              temp1 := temp1*aqnv;
  580.              temp := 2*temp1*root44;
  581.              d4410 := temp*f441*g410;
  582.              d4422 := temp*f442*g422;
  583.              temp1 := temp1*aqnv;
  584.              temp := temp1*root52;
  585.              d5220 := temp*f522*g520;
  586.              d5232 := temp*f523*g532;
  587.              temp := 2*temp1*root54;
  588.              d5421 := temp*f542*g521;
  589.              d5433 := temp*f543*g533;
  590.              xlamo := xmao + xnodeo + xnodeo - thgr - thgr;
  591.              bfact := xlldot + xnodot + xnodot - thdt - thdt;
  592.              bfact := bfact + ssl + ssh + ssh;
  593.              Goto 80;
  594.     { Synchronous resonance terms initialization }
  595.          70: iresfl := 1;
  596.              isynfl := 1;
  597.              g200 := 1 + eqsq*(-2.5 + 0.8125*eqsq);
  598.              g310 := 1 + 2*eqsq;
  599.              g300 := 1 + eqsq*(-6 + 6.60937*eqsq);
  600.              f220 := 0.75*(1 + cosiq)*(1 + cosiq);
  601.              f311 := 0.9375*siniq*siniq*(1 + 3*cosiq) - 0.75*(1 + cosiq);
  602.              f330 := 1 + cosiq;
  603.              f330 := 1.875*f330*f330*f330;
  604.              del1 := 3*xnq*xnq*aqnv*aqnv;
  605.              del2 := 2*del1*f220*g200*q22;
  606.              del3 := 3*del1*f330*g300*q33*aqnv;
  607.              del1 := del1*f311*g310*q31*aqnv;
  608.              fasx2 := 0.13130908;
  609.              fasx4 := 2.8843198;
  610.              fasx6 := 0.37448087;
  611.              xlamo := xmao + xnodeo + omegao - thgr;
  612.              bfact := xlldot + xpidot - thdt;
  613.              bfact := bfact + ssl + ssg + ssh;
  614.          80: xfact := bfact - xnq;
  615.     { Initialize integrator }
  616.              xli := xlamo;
  617.              xni := xnq;
  618.              atime := 0;
  619.              stepp := 720;
  620.              stepn := -720;
  621.              step2 := 259200;
  622.              end; {dpinit}
  623.     dpsec  : begin { Entrance for deep space secular effects }
  624.              xll := xll + ssl*t;
  625.              omgasm := omgasm + ssg*t;
  626.              xnodes := xnodes + ssh*t;
  627.              _em := eo + sse*t;
  628.              xinc := xincl + ssi*t;
  629.              if (xinc >= 0) then
  630.                Goto 90;
  631.              xinc := -xinc;
  632.              xnodes := xnodes  +  pi;
  633.              omgasm := omgasm - pi;
  634.          90: if (iresfl = 0) then
  635.                exit;
  636.         100: if (atime = 0) then
  637.                Goto 170;
  638.              if (t >= 0) and (atime < 0) then
  639.                Goto 170;
  640.              if (t < 0) and (atime >= 0) then
  641.                Goto 170;
  642.         105: if (Abs(t) >= Abs(atime)) then
  643.                Goto 120;
  644.              delt := stepp;
  645.              if (t >= 0) then
  646.                delt := stepn;
  647.         110: iret := 100; {assign 100 to iret}
  648.              Goto 160;
  649.         120: delt := stepn;
  650.              if (t > 0) then
  651.                delt := stepp;
  652.         125: if (Abs(t - atime) < stepp) then
  653.                Goto 130;
  654.              iret := 125; {assign 125 to iret}
  655.              Goto 160;
  656.         130: ft := t - atime;
  657.              iretn := 140; {assign 140 to iretn}
  658.              Goto 150;
  659.         140: xn := xni + xndot*ft + xnddt*ft*ft*0.5;
  660.              xl := xli + xldot*ft + xndot*ft*ft*0.5;
  661.              temp := -xnodes + thgr + t*thdt;
  662.              xll := xl - omgasm + temp;
  663.              if (isynfl = 0) then
  664.                xll := xl + temp + temp;
  665.              exit;
  666.     { Dot terms calculated }
  667.         150: if (isynfl = 0) then
  668.                Goto 152;
  669.              xndot := del1*Sin(xli - fasx2) + del2*Sin(2*(xli - fasx4))
  670.                     + del3*Sin(3*(xli - fasx6));
  671.              xnddt := del1*Cos(xli - fasx2)
  672.                     + 2*del2*Cos(2*(xli - fasx4))
  673.                     + 3*del3*Cos(3*(xli - fasx6));
  674.              Goto 154;
  675.         152: xomi := omegaq + omgdt*atime;
  676.              x2omi := xomi + xomi;
  677.              x2li := xli + xli;
  678.              xndot := d2201*Sin(x2omi + xli - g22)
  679.                     + d2211*Sin(xli - g22)
  680.                     + d3210*Sin(xomi + xli - g32)
  681.                     + d3222*Sin(-xomi + xli - g32)
  682.                     + d4410*Sin(x2omi + x2li - g44)
  683.                     + d4422*Sin(x2li - g44)
  684.                     + d5220*Sin(xomi + xli - g52)
  685.                     + d5232*Sin(-xomi + xli - g52)
  686.                     + d5421*Sin(xomi + x2li - g54)
  687.                     + d5433*Sin(-xomi + x2li - g54);
  688.              xnddt := d2201*Cos(x2omi + xli - g22)
  689.                     + d2211*Cos(xli - g22)
  690.                     + d3210*Cos(xomi + xli - g32)
  691.                     + d3222*Cos(-xomi + xli - g32)
  692.                     + d5220*Cos(xomi + xli - g52)
  693.                     + d5232*Cos(-xomi + xli - g52)
  694.                     + 2*(d4410*Cos(x2omi + x2li - g44)
  695.                     + d4422*Cos(x2li - g44)
  696.                     + d5421*Cos(xomi + x2li - g54)
  697.                     + d5433*Cos(-xomi + x2li - g54));
  698.         154: xldot := xni + xfact;
  699.              xnddt := xnddt*xldot;
  700.              case iretn of
  701.                140 : Goto 140;
  702.                165 : Goto 165;
  703.              else
  704.                      Halt;
  705.              end; {case}
  706.     { Integrator }
  707.         160: iretn := 165; {assign 165 to iretn}
  708.              Goto 150;
  709.         165: xli := xli + xldot*delt + xndot*step2;
  710.              xni := xni + xndot*delt + xnddt*step2;
  711.              atime := atime + delt;
  712.              case iret of
  713.                100 : Goto 100;
  714.                125 : Goto 125;
  715.              else
  716.                      Halt;
  717.              end; {case}
  718.     { Epoch restart }
  719.         170: if (t >= 0) then
  720.                Goto 175;
  721.              delt := stepn;
  722.              Goto 180;
  723.         175: delt := stepp;
  724.         180: atime := 0;
  725.              xni := xnq;
  726.              xli := xlamo;
  727.              Goto 125;
  728.              end; {dpsec}
  729.     dpper  : begin { Entrance for lunar-solar periodics }
  730.              sinis := Sin(xinc);
  731.              cosis := Cos(xinc);
  732.              if (Abs(savtsn - t) < 30) then
  733.                Goto 210;
  734.              savtsn := t;
  735.              zm := zmos + zns*t;
  736.         205: zf := zm + 2*zes*Sin(zm);
  737.              sinzf := Sin(zf);
  738.              f2 := 0.5*sinzf*sinzf - 0.25;
  739.              f3 := -0.5*sinzf*Cos(zf);
  740.              ses := se2*f2 + se3*f3;
  741.              sis := si2*f2 + si3*f3;
  742.              sls := sl2*f2 + sl3*f3 + sl4*sinzf;
  743.              sghs := sgh2*f2 + sgh3*f3 + sgh4*sinzf;
  744.              shs := sh2*f2 + sh3*f3;
  745.              zm := zmol + znl*t;
  746.              zf := zm + 2*zel*Sin(zm);
  747.              sinzf := Sin(zf);
  748.              f2 := 0.5*sinzf*sinzf - 0.25;
  749.              f3 := -0.5*sinzf*Cos(zf);
  750.              sel := ee2*f2 + e3*f3;
  751.              sil := xi2*f2 + xi3*f3;
  752.              sll := xl2*f2 + xl3*f3 + xl4*sinzf;
  753.              sghl := xgh2*f2 + xgh3*f3 + xgh4*sinzf;
  754.              sh1 := xh2*f2 + xh3*f3;
  755.              pe := ses + sel;
  756.              pinc := sis + sil;
  757.              pl := sls + sll;
  758.         210: pgh := sghs + sghl;
  759.              ph := shs + sh1;
  760.              xinc := xinc + pinc;
  761.              _em := _em + pe;
  762.              if (xqncl < 0.2) then
  763.                Goto 220;
  764.              Goto 218;
  765.     { Apply periodics directly }
  766.         218: ph := ph/siniq;
  767.              pgh := pgh - cosiq*ph;
  768.              omgasm := omgasm + pgh;
  769.              xnodes := xnodes + ph;
  770.              xll := xll + pl;
  771.              Goto 230;
  772.     { Apply periodics with Lyddane modification }
  773.         220: sinok := Sin(xnodes);
  774.              cosok := Cos(xnodes);
  775.              alfdp := sinis*sinok;
  776.              betdp := sinis*cosok;
  777.              dalf := ph*cosok + pinc*cosis*sinok;
  778.              dbet := -ph*sinok + pinc*cosis*cosok;
  779.              alfdp := alfdp + dalf;
  780.              betdp := betdp + dbet;
  781.              xls := xll + omgasm + cosis*xnodes;
  782.              dls := pl + pgh - pinc*xnodes*sinis;
  783.              xls := xls + dls;
  784.              xnodes := Actan(alfdp,betdp);
  785.              xll := xll + pl;
  786.              omgasm := xls - xll - Cos(xinc)*xnodes;
  787.         230: end; {dpper}
  788.     end; {case}
  789.   end; {Procedure Deep}
  790.  
  791. Procedure Call_dpinit(var eosq,sinio,cosio,betao,aodp,theta2,
  792.                           sing,cosg,betao2,xmdot,omgdot,
  793.                           xnodott,xnodpp : double);
  794.   begin
  795.   eqsq   := eosq;    siniq  := sinio;   cosiq  := cosio;   rteqsq := betao;
  796.   ao     := aodp;    cosq2  := theta2;  sinomo := sing;    cosomo := cosg;
  797.   bsq    := betao2;  xlldot := xmdot;   omgdt  := omgdot;  xnodot := xnodott;
  798.   xnodp  := xnodpp;
  799.   Deep(1);
  800.   eosq   := eqsq;    sinio  := siniq;   cosio  := cosiq;   betao  := rteqsq;
  801.   aodp   := ao;      theta2 := cosq2;   sing   := sinomo;  cosg   := cosomo;
  802.   betao2 := bsq;     xmdot  := xlldot;  omgdot := omgdt;   xnodott := xnodot;
  803.   xnodpp  := xnodp;
  804.   end; {Procedure Call_dpinit}
  805.  
  806. Procedure Call_dpsec(var xmdf,omgadf,xnode,emm,xincc,xnn,tsince : double);
  807.   begin
  808.   xll    := xmdf;    omgasm := omgadf;  xnodes := xnode;   {_em     := emm;
  809.   xinc   := xincc;}   xn     := xnn;     t      := tsince;
  810.   Deep(2);
  811.   xmdf   := xll;     omgadf := omgasm;  xnode  := xnodes;  emm    := _em;
  812.   xincc  := xinc;    xnn    := xn;      tsince := t;
  813.   end; {Procedure Call_dpsec}
  814.  
  815. Procedure Call_dpper(var e,xincc,omgadf,xnode,xmam : double);
  816.   begin
  817.   _em     := e;       xinc   := xincc;   omgasm := omgadf;  xnodes := xnode;
  818.   xll    := xmam;
  819.   Deep(3);
  820.   e      := _em;      xincc  := xinc;    omgadf := omgasm;  xnode  := xnodes;
  821.   xmam   := xll;
  822.   end; {Procedure Call_dpper}
  823.  
  824. Procedure SDP4(tsince : double;
  825.             var iflag : integer;
  826.           var pos,vel : vector);
  827.   label
  828.     9,10,90,100,130,140;
  829.   const
  830.     a1     : double = 0;  a3ovk2 : double = 0;  ao     : double = 0;
  831.     aodp   : double = 0;  aycof  : double = 0;  betao  : double = 0;
  832.     betao2 : double = 0;  c1     : double = 0;  c2     : double = 0;
  833.     c4     : double = 0;  coef   : double = 0;  coef1  : double = 0;
  834.     cosg   : double = 0;  cosio  : double = 0;  del1   : double = 0;
  835.     delo   : double = 0;  eeta   : double = 0;  eosq   : double = 0;
  836.     eta    : double = 0;  etasq  : double = 0;  omgdot : double = 0;
  837.     perige : double = 0;  pinvsq : double = 0;  psisq  : double = 0;
  838.     qoms24 : double = 0;  s4     : double = 0;  sing   : double = 0;
  839.     sinio  : double = 0;  t2cof  : double = 0;  temp1  : double = 0;
  840.     temp2  : double = 0;  temp3  : double = 0;  theta2 : double = 0;
  841.     theta4 : double = 0;  tsi    : double = 0;  x1m5th : double = 0;
  842.     x1mth2 : double = 0;  x3thm1 : double = 0;  x7thm1 : double = 0;
  843.     xhdot1 : double = 0;  xlcof  : double = 0;  xmdot  : double = 0;
  844.     xnodcf : double = 0;  xnodot : double = 0;  xnodp  : double = 0;
  845.   var
  846.     i                                                  : integer;
  847.     a,axn,ayn,aynl,beta,betal,capu,cos2u,cosepw,cosik,
  848.     cosnok,cosu,cosuk,e,ecose,elsq,em,epw,esine,omgadf,
  849.     pl,r,rdot,rdotk,rfdot,rfdotk,rk,sin2u,sinepw,sinik,
  850.     sinnok,sinu,sinuk,temp,temp4,temp5,temp6,tempa,
  851.     tempe,templ,tsq,u,uk,ux,uy,uz,vx,vy,vz,xinc,xinck,
  852.     xl,xll,xlt,xmam,xmdf,xmx,xmy,xn,xnoddf,xnode,xnodek,
  853.     x,y,z,xdot,ydot,zdot                               : double;
  854.   begin
  855.   if (iflag = 0) then
  856.     goto 100;
  857. { Recover original mean motion (xnodp) and semimajor axis (aodp) }
  858. { from input elements. }
  859.   a1 := Power(xke/xno,tothrd);
  860.   cosio := cos(xincl);
  861.   theta2 := cosio*cosio;
  862.   x3thm1 := 3*theta2 - 1;
  863.   eosq := eo*eo;
  864.   betao2 := 1 - eosq;
  865.   betao := sqrt(betao2);
  866.   del1 := 1.5*ck2*x3thm1/(a1*a1*betao*betao2);
  867.   ao := a1*(1 - del1*(0.5*tothrd + del1*(1 + 134/81*del1)));
  868.   delo := 1.5*ck2*x3thm1/(ao*ao*betao*betao2);
  869.   xnodp := xno/(1 + delo);
  870.   aodp := ao/(1 - delo);
  871. { Initialization }
  872. { For perigee below 156 km, the values of s and qoms2t are altered. }
  873.   s4 := s;
  874.   qoms24 := qoms2t;
  875.   perige := (aodp*(1 - eo) - ae)*xkmper;
  876.   if (perige >= 156) then goto 10;
  877.   s4 := perige - 78;
  878.   if (perige > 98) then goto 9;
  879.   s4 := 20;
  880.   9:
  881.   qoms24 := Power((120 - s4)*ae/xkmper,4);
  882.   s4 := s4/xkmper + ae;
  883.  10:
  884.   pinvsq := 1/(aodp*aodp*betao2*betao2);
  885.   sing := sin(omegao);
  886.   cosg := cos(omegao);
  887.   tsi := 1/(aodp - s4);
  888.   eta := aodp*eo*tsi;
  889.   etasq := eta*eta;
  890.   eeta := eo*eta;
  891.   psisq := abs(1 - etasq);
  892.   coef := qoms24*Power(tsi,4);
  893.   coef1 := coef/Power(psisq,3.5);
  894.   c2 := coef1*xnodp*(aodp*(1 + 1.5*etasq + eeta*(4 + etasq))
  895.       + 0.75*ck2*tsi/psisq*x3thm1*(8 + 3*etasq*(8 + etasq)));
  896.   c1 := bstar*c2;
  897.   sinio := sin(xincl);
  898.   a3ovk2 := -xj3/ck2*Power(ae,3);
  899.   x1mth2 := 1 - theta2;
  900.   c4 := 2*xnodp*coef1*aodp*betao2*(eta*(2 + 0.5*etasq)
  901.       + eo*(0.5 + 2*etasq) - 2*ck2*tsi/(aodp*psisq)
  902.       *(-3*x3thm1*(1 - 2*eeta + etasq*(1.5 - 0.5*eeta))
  903.       + 0.75*x1mth2*(2*etasq - eeta*(1 + etasq))*cos(2*omegao)));
  904.   theta4 := theta2*theta2;
  905.   temp1 := 3*ck2*pinvsq*xnodp;
  906.   temp2 := temp1*ck2*pinvsq;
  907.   temp3 := 1.25*ck4*pinvsq*pinvsq*xnodp;
  908.   xmdot := xnodp + 0.5*temp1*betao*x3thm1
  909.          + 0.0625*temp2*betao*(13 - 78*theta2 + 137*theta4);
  910.   x1m5th := 1 - 5*theta2;
  911.   omgdot := -0.5*temp1*x1m5th + 0.0625*temp2*(7 - 114*theta2 + 395*theta4)
  912.           + temp3*(3 - 36*theta2 + 49*theta4);
  913.   xhdot1 := -temp1*cosio;
  914.   xnodot := xhdot1 + (0.5*temp2*(4 - 19*theta2)
  915.           + 2*temp3*(3 - 7*theta2))*cosio;
  916.   xnodcf := 3.5*betao2*xhdot1*c1;
  917.   t2cof := 1.5*c1;
  918.   xlcof := 0.125*a3ovk2*sinio*(3 + 5*cosio)/(1 + cosio);
  919.   aycof := 0.25*a3ovk2*sinio;
  920.   x7thm1 := 7*theta2 - 1;
  921.  90:
  922.   iflag := 0;
  923.   Call_dpinit(eosq,sinio,cosio,betao,aodp,theta2,sing,cosg,
  924.               betao2,xmdot,omgdot,xnodot,xnodp);
  925. { Update for secular gravity and atmospheric drag }
  926. 100:
  927.   xmdf := xmo + xmdot*tsince;
  928.   omgadf := omegao + omgdot*tsince;
  929.   xnoddf := xnodeo + xnodot*tsince;
  930.   tsq := tsince*tsince;
  931.   xnode := xnoddf + xnodcf*tsq;
  932.   tempa := 1 - c1*tsince;
  933.   tempe := bstar*c4*tsince;
  934.   templ := t2cof*tsq;
  935.   xn := xnodp;
  936.   Call_dpsec(xmdf,omgadf,xnode,em,xinc,xn,tsince);
  937.   a := Power(xke/xn,tothrd)*Sqr(tempa);
  938.   e := em - tempe;
  939.   xmam := xmdf + xnodp*templ;
  940.   Call_dpper(e,xinc,omgadf,xnode,xmam);
  941.   xl := xmam + omgadf + xnode;
  942.   beta := sqrt(1 - e*e);
  943.   xn := xke/Power(a,1.5);
  944. { Long period periodics }
  945.   axn := e*cos(omgadf);
  946.   temp := 1/(a*beta*beta);
  947.   xll := temp*xlcof*axn;
  948.   aynl := temp*aycof;
  949.   xlt := xl + xll;
  950.   ayn := e*sin(omgadf) + aynl;
  951. { Solve Kepler's Equation }
  952.   capu := fmod2p(xlt - xnode);
  953.   temp2 := capu;
  954.   for i := 1 to 10 do
  955.     begin
  956.     sinepw := sin(temp2);
  957.     cosepw := cos(temp2);
  958.     temp3 := axn*sinepw;
  959.     temp4 := ayn*cosepw;
  960.     temp5 := axn*cosepw;
  961.     temp6 := ayn*sinepw;
  962.     epw := (capu - temp4 + temp3 - temp2)/(1 - temp5 - temp6) + temp2;
  963.     if (abs(epw - temp2) <= e6a) then goto 140;
  964. 130:
  965.     temp2 := epw;
  966.     end; {for i}
  967. { Short period preliminary quantities }
  968. 140:
  969.   ecose := temp5 + temp6;
  970.   esine := temp3 - temp4;
  971.   elsq := axn*axn + ayn*ayn;
  972.   temp := 1 - elsq;
  973.   pl := a*temp;
  974.   r := a*(1 - ecose);
  975.   temp1 := 1/r;
  976.   rdot := xke*sqrt(a)*esine*temp1;
  977.   rfdot := xke*sqrt(pl)*temp1;
  978.   temp2 := a*temp1;
  979.   betal := sqrt(temp);
  980.   temp3 := 1/(1 + betal);
  981.   cosu := temp2*(cosepw - axn + ayn*esine*temp3);
  982.   sinu := temp2*(sinepw - ayn - axn*esine*temp3);
  983.   u := actan(sinu,cosu);
  984.   sin2u := 2*sinu*cosu;
  985.   cos2u := 2*cosu*cosu - 1;
  986.   temp := 1/pl;
  987.   temp1 := ck2*temp;
  988.   temp2 := temp1*temp;
  989. { Update for short periodics }
  990.   rk := r*(1 - 1.5*temp2*betal*x3thm1) + 0.5*temp1*x1mth2*cos2u;
  991.   uk := u - 0.25*temp2*x7thm1*sin2u;
  992.   xnodek := xnode + 1.5*temp2*cosio*sin2u;
  993.   xinck := xinc + 1.5*temp2*cosio*sinio*cos2u;
  994.   rdotk := rdot - xn*temp1*x1mth2*sin2u;
  995.   rfdotk := rfdot + xn*temp1*(x1mth2*cos2u + 1.5*x3thm1);
  996. { Orientation vectors }
  997.   sinuk := sin(uk);
  998.   cosuk := cos(uk);
  999.   sinik := sin(xinck);
  1000.   cosik := cos(xinck);
  1001.   sinnok := sin(xnodek);
  1002.   cosnok := cos(xnodek);
  1003.   xmx := -sinnok*cosik;
  1004.   xmy := cosnok*cosik;
  1005.   ux := xmx*sinuk + cosnok*cosuk;
  1006.   uy := xmy*sinuk + sinnok*cosuk;
  1007.   uz := sinik*sinuk;
  1008.   vx := xmx*cosuk - cosnok*sinuk;
  1009.   vy := xmy*cosuk - sinnok*sinuk;
  1010.   vz := sinik*cosuk;
  1011. { Position and velocity }
  1012.   x := rk*ux;  pos[1] := x;
  1013.   y := rk*uy;  pos[2] := y;
  1014.   z := rk*uz;  pos[3] := z;
  1015.   xdot := rdotk*ux + rfdotk*vx;  vel[1] := xdot;
  1016.   ydot := rdotk*uy + rfdotk*vy;  vel[2] := ydot;
  1017.   zdot := rdotk*uz + rfdotk*vz;  vel[3] := zdot;
  1018.   end; {Procedure SDP4}
  1019.  
  1020. Procedure SGP(time : double;
  1021.        var pos,vel : vector);
  1022.   var
  1023.     tsince : double;
  1024.   begin
  1025.   tsince := (time - julian_epoch) * xmnpda;
  1026.   if ideep = 0 then
  1027.     SGP4(tsince,iflag,pos,vel)
  1028.   else
  1029.     SDP4(tsince,iflag,pos,vel);
  1030.   end; {Procedure SGP}
  1031.  
  1032.   begin
  1033.   Define_Derived_Constants;
  1034.  
  1035. end.
  1036.